library(foreign)
library(ggplot2)
library(grid)
library(gridExtra)
library(scales)

rm(list=ls())

setwd("Yourpath")

source("fine_grid.R") # ggplot layers
source("summarySE.R") # summary stats (for figure 5)


########################################################################
#                                                                      #
#      PLOT MAIN RESULTS SYNTHETIC CONTROL METHOD (Figure 1)           #
#                                                                      #
########################################################################



# (i) Levels (Figure 1a)

Y_0 = read.csv("Y_synthetic_nonannual.txt", sep = "\t")[,2]*100
Y_1 = read.csv("Y_treated_nonannual.txt", sep = "\t")[,2]*100
X=read.csv("Y_treated_nonannual.txt", sep = "\t")[,1] 
W_weights_nonannual=read.csv("W_weights_nonannual.txt", sep = "\t")[,3] 


yearcutoffs <- c(1,11,21,42,59,74,106,127,130)

X_equalspace <- matrix(NA,length(X),1)

for (i in 2:length(yearcutoffs)){
X_min <- min(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
X_max <- max(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
X_dist <- X_max-X_min 
X_equalspace[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]  <- (i-1)+ (X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]-X_min)/X_dist 
}

X_equal_10years <- X_equalspace[is.element(X,yearcutoffs)]
X_equalspace[c(128:130)] <- (length(yearcutoffs)-1)+(X_equalspace[127]-X_equalspace[126])*c(1:3)



gray.scale = .7
gray.scale2 = .7
dim.factor = .1
Legend = c("Vaud                  ","Synthetic Vaud")
Ylab = c("Turnout (%)")
Xlab = c("Year") 
Ylim = c(0,100) 
Xlim = c(min(X_equalspace),max(X_equalspace))
gap <- Y_1 - Y_0
plot(X_equalspace, Y_1, t = "l", col = "white", 
     lwd = 1,  ylab = Ylab, xlab = Xlab, ylim = Ylim, xlim=Xlim,
     yaxs = "i",xaxt="n") 
polygon(x=c(3.571429,3.571429,4.941176,4.941176),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE)
polygon(x=c( 5.466667, 5.466667,5.733333,5.733333),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE) 
lines(X_equalspace, Y_1, lwd = 1)
lines(X_equalspace, Y_0, col = gray(.55), lwd = 1)
points(X_equalspace, Y_1, pch=19, cex = 4/7 )
points(X_equalspace, Y_0, pch=19, col = gray(.55), cex = 4/7)

axis(1,las=1,labels=c(seq(1900,1970,10)),at=X_equal_10years[1:8])

legend(c("topright"), legend = Legend, lty = 1:1, col = c("black", 
                                                          gray(.55)), lwd = c(1.75, 1.75), cex = 6/7)

dev.copy2pdf(file="Figure1a.pdf",width=8,height=6)


# (ii) Effect (Figure 1b)

Y_0 = read.csv("Y_synthetic_nonannual.txt", sep = "\t")[,2]*100
Y_1 = read.csv("Y_treated_nonannual.txt", sep = "\t")[,2]*100
X=read.csv("Y_treated_nonannual.txt", sep = "\t")[,1] 
W_weights=read.csv("W_weights_nonannual.txt", sep = "\t")[,3] 



gap =Y_1-Y_0
gray.scale = .7
gray.scale2 = .7
dim.factor = .3
Ylab = c("Effect on Turnout (Percentage Points)")
Xlab = c("Year") 
Ylim = c(-17,65) 
Xlim = c(min(X_equalspace),max(X_equalspace))
gap <- Y_1 - Y_0
plot(X_equalspace, gap, t = "l", col = "white", 
     lwd = 1,  ylab = Ylab, xlab = Xlab, ylim = Ylim, xlim=Xlim,
     yaxs = "i",xaxt="n") 
polygon(x=c(3.571429,3.571429,4.941176,4.941176),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE)
polygon(x=c( 5.466667, 5.466667,5.733333,5.733333),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE) 
lines(X_equalspace, gap, lwd = 1)
points(X_equalspace, gap, pch=19, cex = 4/7 )

axis(1,las=1,labels=c(seq(1900,1970,10)),at=X_equal_10years[1:8])

dev.copy2pdf(file="Figure1b.pdf",width=8,height=6)




################################################
#                                              #
#       PLACEBO TEST IN TIME (Figure 4a)       #
#                                              #
################################################



Y_0 = read.csv("Y_synthetic_placebotime_nonannual.txt", sep = "\t")[,2]*100
Y_1 = read.csv("Y_treated_placebotime_nonannual.txt", sep = "\t")[,2]*100
X=read.csv("Y_treated_nonannual.txt", sep = "\t")[,1]


yearcutoffs <- c(1,11,21,42,59,74,106,127,130)

X_equalspace <- matrix(NA,length(X),1)

for (i in 2:length(yearcutoffs)){
  X_min <- min(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
  X_max <- max(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
  X_dist <- X_max-X_min 
  X_equalspace[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]  <- (i-1)+ (X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]-X_min)/X_dist 
}

X_equal_10years <- X_equalspace[is.element(X,yearcutoffs)]

X_equalspace[c(128:130)] <- (length(yearcutoffs)-1)+(X_equalspace[127]-X_equalspace[126])*c(1:3)

Y_0 = read.csv("Y_synthetic_placebotime_nonannual.txt", sep = "\t")[,2]*100
Y_1 = read.csv("Y_treated_placebotime_nonannual.txt", sep = "\t")[,2]*100
X=read.csv("Y_treated_nonannual.txt", sep = "\t")[,1] 


yearcutoffs <- c(1,11,21,42,59,74,106,127,130)
X_equalspace <- matrix(NA,length(X),1)
for (i in 2:length(yearcutoffs)){
  X_min <- min(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
  X_max <- max(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
  X_dist <- X_max-X_min 
  X_equalspace[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]  <- (i-1)+ (X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]-X_min)/X_dist 
}

X_equal_10years <- X_equalspace[is.element(X,yearcutoffs)]

X_equalspace[c(128:130)] <- (length(yearcutoffs)-1)+(X_equalspace[127]-X_equalspace[126])*c(1:3)


gap =Y_1-Y_0


gray.scale = .7
gray.scale2 = .9
dim.factor = .3
Ylab = c("Effect on Turnout (Percentage Points)")
Xlab = c("Year") 
Ylim = c(-40,72) 
gap <- Y_1 - Y_0
plot(X_equalspace, gap, t = "l", col = "white", 
     lwd = 1, ylab = Ylab, xlab = Xlab, ylim = Ylim, 
     xaxt = "n", yaxs = "i") 
polygon(x=c(2.4,2.4,3.571429,3.571429),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale2),border=gray(gray.scale2),lty=1, bg=TRUE)
polygon(x=c(3.571429,3.571429,4.941176,4.941176),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE)
polygon(x=c( 5.466667, 5.466667,5.733333,5.733333),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE) 
abline(h=0, lwd=1, lty = "dashed")      

lines(X_equalspace, gap, lwd = 1)
points(X_equalspace, gap, pch=19, cex = 4/7)


axis(1,las=1,labels=c(seq(1900,1970,10)),at=X_equal_10years[1:8])

dev.copy2pdf(file="Figure4a.pdf",width=8,height=6)



################################################
#                                              #
#      PLACEBO TEST IN SPACE  (FIGURE 4B)      #
#                                              #
################################################


Y_0 = read.csv("Y_synthetic_nonannual.txt", sep = "\t")[,2]*100
Y_1 = read.csv("Y_treated_nonannual.txt", sep = "\t")[,2]  *100
X=read.csv("Y_treated_nonannual.txt", sep = "\t")[,1] 
gray.scale = .7
dim.factor = .3
Ylab = c("Effect on Turnout (Percentage Points)")
Xlab = c("Year") 
Ylim = c(-40,72) 

fs <- c(1,11,21,42,59,74,106,127,130)
X_equalspace <- matrix(NA,length(X),1)
for (i in 2:length(yearcutoffs)){
  X_min <- min(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
  X_max <- max(X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])])
  X_dist <- X_max-X_min 
  X_equalspace[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]  <- (i-1)+ (X[is.element(X,yearcutoffs[i]:yearcutoffs[(i-1)])]-X_min)/X_dist 
}

X_equal_10years <- X_equalspace[is.element(X,yearcutoffs)]

X_equalspace <- X_equalspace[c(1:127)] # cut off 1970 referendums
X_equalspace[c(128:130)] <- (length(yearcutoffs)-1)+(X_equalspace[127]-X_equalspace[126])*c(1:3)


gap =Y_1-Y_0

plot(X_equalspace, gap, t = "l", col = "white", 
     lwd = 1, ylab = Ylab, xlab = Xlab, ylim = Ylim, 
     xaxs = "i", yaxs = "i", xaxt = "n") 
polygon(x=c(3.571429,3.571429,4.941176,4.941176),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE)
polygon(x=c( 5.466667, 5.466667,5.733333,5.733333),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE) 
abline(h=0, lwd=1, lty = "dashed")   

for (i in c(1:14)) {
  Y_0 = read.csv(paste("Y_synthetic_nonannual_",i,".txt", sep = ""), sep = "\t")[,2]*100
  Y_1 = read.csv(paste("Y_treated_nonannual_",i,".txt", sep = ""), sep = "\t")[,2]*100
  X_equalspace <- X_equalspace[c(1:127)] # cut off 1970 referendums
  Y_0 <- Y_0[c(1:127)]
  Y_1 <- Y_1[c(1:127)] 
  gap <- Y_1 - Y_0  
  lines(X_equalspace, gap, pch=19, cex = 4/7, col = "grey" )
}


Y_0 = read.csv("Y_synthetic_nonannual.txt", sep = "\t")[,2]
Y_1 = read.csv("Y_treated_nonannual.txt", sep = "\t")[,2]  
Y_0 <- Y_0[c(1:127)]*100
Y_1 <- Y_1[c(1:127)]*100
gap <- Y_1 - Y_0  

lines(X_equalspace, gap, pch=19, cex = 4/7, col = "black" )
points(X_equalspace, gap, pch=19, cex = 4/7)

axis(1,las=1,labels=c(seq(1900,1970,10)),at=X_equal_10years[1:8])


dev.copy2pdf(file="Figure4b.pdf",width=8,height=6)



########################################################################
#                                                                      #
#      PLOT USING ANNUAL SYNTHETIC CONTROL METHOD WEIGHTS (Figure 7)   #
#                                                                      #
########################################################################



mean(gap[X_equalspace>3.57142& X_equalspace<4.94117| X_equalspace>5.466 & X_equalspace<5.73])

# compute ATET using weights for annual weights 


W_weights=read.csv("W_weights.txt", sep = "\t")[,3]  # annual weights
data_refs <- read.dta("federal_referendums_1900_1970_covariates_wide.dta") # non-annual data

data_refs_vd <- data_refs[data_refs$id==22,] # data for VD
data_refs <- data_refs[data_refs$id!=22,] # exclude VD
data_refs <- data_refs[,c(2:dim(data_refs)[2])] # exclude VD
data_refs_vd <- data_refs_vd[,c(2:dim(data_refs_vd)[2])] # exclude VD

gap_alt_weights <- data_refs_vd-as.vector(W_weights)%*%as.matrix(data_refs)


mean(gap_alt_weights[X_equalspace>3.57142& X_equalspace<4.94117| X_equalspace>5.466 & X_equalspace<5.73])*100


# compute ATET using weights for non-annual weights 


Y_0 = read.csv("Y_synthetic_nonannual.txt", sep = "\t")[,2]*100
Y_1 = read.csv("Y_treated_nonannual.txt", sep = "\t")[,2]*100
Y_0_alt=as.vector(W_weights)%*%as.matrix(data_refs)*100
X=read.csv("Y_treated_nonannual.txt", sep = "\t")[,1] 


X_equal_10years <- X_equalspace[is.element(X,yearcutoffs)]
X_equalspace[c(128:130)] <- (length(yearcutoffs)-1)+(X_equalspace[127]-X_equalspace[126])*c(1:3)


gap =Y_1-Y_0
gray.scale = .7
gray.scale2 = .7
dim.factor = .3
Legend = c("Vaud                  ","Synth. Vaud (RDT)","Synth. Vaud (ATA)")

Ylab = c("Referendum Day Turnout (%)")
Xlab = c("Year") 
Ylim = c(0,100) 
Xlim = c(min(X_equalspace),max(X_equalspace))
gap <- Y_1 - Y_0
plot(X_equalspace, gap, t = "l", col = "white", 
     lwd = 1,  ylab = Ylab, xlab = Xlab, ylim = Ylim, xlim=Xlim,
     yaxs = "i",xaxt="n") 
polygon(x=c(3.571429,3.571429,4.941176,4.941176),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE)
polygon(x=c( 5.466667, 5.466667,5.733333,5.733333),
        y=c(Ylim[1]+dim.factor,Ylim[2]-dim.factor,Ylim[2]-dim.factor,Ylim[1]+dim.factor),
        col=gray(gray.scale),border=gray(gray.scale),lty=1, bg=TRUE) 

lines(X_equalspace, Y_1, lwd = 1)
lines(X_equalspace, as.numeric(Y_0_alt), col = gray(.55), lwd = 1)
lines(X_equalspace, as.numeric(Y_0), col = gray(.55), lwd = 1)

points(X_equalspace, Y_1, pch=19, cex =5/7 )
points(X_equalspace, as.numeric(Y_0_alt), pch=5, col = gray(.55), cex = 5/7)
points(X_equalspace, as.numeric(Y_0), pch=19, col = gray(.55), cex = 5/7)

axis(1,las=1,labels=c(seq(1900,1970,10)),at=X_equal_10years[1:8])

legend(c("topright"), legend = Legend, lty = 1:1, pch=c(19,19,5),col = c("black", 
                                                                         gray(.55), gray(.55)), lwd = c(1.75, 1.75,1.75), cex = 5/7)

axis(1,las=1,labels=c(seq(1900,1970,10)),at=X_equal_10years[1:8])
dev.copy2pdf(file="Figure7.pdf",width=8,height=6)
